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Over  generations,  a  species  will  evolve  to  increase  its  fitness  to  its  environment.  Natural 
selection  is  the  primary  mechanism  by  which  this  occurs,  wherein  environmental  pressures 
cause  individuals  with  more  fit  phenotypes  to  reproduce  at  a  higher  rate.  For  the  purpose 
of  modeling,  the  space  of  possible  heritable  traits  may  be  viewed  as  a  continuum,  and  a 
particular  individual  as  inhabiting  a  point  in  trait  space.  Making  the  biomass  continuous  as 
well,  a  population  is  viewed  as  a  distribution  in  trait  space.  In  this  view,  natural  selection  is 
the  process  by  which  the  distribution  of  biomass  propagates  through  trait  space,  generally 
increasing  its  average  fitness. 

1  Building  blocks  for  models 

The  first  step  in  modeling  a  system  in  which  the  dominant  traits  can  change  is  defining  a 
suitable  measure  of  how  much  of  a  species  there  is.  While  the  actual  number  of  members  of 
the  species  is  an  option,  a  more  useful  measure  is  the  biomass,  b,  which  readily  translates 
across  species.  The  biomass  density  can  be  treated  as  a  function  over  space,  time,  and  this 
“trait-space” 

fr(x,  t\s\,  S2,  •••) 

where  the  variables  s  =  {si,  S2,  •  •  •}  specify  the  phenotype,  or  traits,  of  the  organism.  In 
addition  to  the  biomass,  the  rates  of  survival  and  reproduction  depend  on  the  environmental 
conditions  and  vary  with  phenotype.  These  rates  are  defined  as 

n  =  n(s\E)  ,  n{Sl\E)  y  tz(s2\e) 

where  E  represents  the  environmental  conditions.  The  biomass  evolution  is  related  to  the 
growth  rate  (including  reproduction)  by 

|&(s)  =  ft(s|£)6(s)...  (1) 

indicating  that  new  organisms  (represented  as  biomass)  have  the  same  values  of  s  as  their 
parents.  However,  in  nature  genetic  mutation  can  occur  causing  the  offspring  of  parents  to 
have  a  slightly  different  trait  values.  (And,  since  these  represent  phenotypes,  there  is  some 
degree  of  natural  variation  as  well.)  This  can  be  modeled  by  changing  Equation  1  to 

|U(s)  =  K(s\E)b(s)  +  V2smb  (2) 
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where  m  is  a  type  of  trait  diffusion  (which  can  be  proportional  to  the  reproduction  rate). 
The  rate  appears  within  the  derivatives  for  the  same  reasons  it  does  in  the  Fokker-Planck 
equation  (or  the  representation  of  kinesis  in  a  previous  lecture)  -  the  length  of  the  random 
jump  in  trait  space  is  set  at  the  beginning  of  the  excursion  (and  there  is  no  reason  to  expect 
non-divergence).  In  any  case,  we  shall  just  use  the  over-simplification  of  treating  m  as 
constant. 

The  simplified  model  created  in  Equation  2  obviously  has  a  number  of  flaws.  If  a  trait 
is  not  represented  in  the  initial  modeling  it  will  never  arise  long  term  meaning  that  a  full 
trait  space  needs  to  be  known  a  priori.  Additionally,  the  time  scales  are  not  well-defined. 
While  there  is  some  idea  of  the  rate  of  genetic  changes,  not  much  about  how  this  translates 
into  alterations  of  the  phenotype  is  understood.  Finally,  since  natural  selection  acts  on  an 
individual  basis;  a  model  based  on  biomass  may  not  accurately  model  the  dynamics  for  low 
values  where  extinction  is  probable. 

If  we  extend  the  above  model  to  start  accounting  for  advection  and  diffusion  in  space  as 
well  as  diffusion  in  trait  space,  we  can  use  the  following  equation  to  model  a  wide  variety 
of  systems 

-  VuVb  =  bK(s\E)  +  V2mb. 

With  this  governing  equation  the  following  topics  will  be  examined 

•  Evolutionarily  stable  strategies:  what  s  values  are  selected  for? 

•  Relationship  to  stability  theory 

•  Adaptive  dynamics  in  the  presence  of  mutation 

•  Effects  of  mixing,  diffusion,  and  advection 


2  A  Simple  Model 


To  make  a  model  which  can  be  studied  easily,  it  is  necessary  to  create  a  function  for 
the  growth  rate.  For  a  first  model,  we  make  a  number  of  simplifying  assumptions:  the 
organisms  will  have  a  ID  trait  space,  s  =  s,  and  use  a  common  resource,  N,  so  that  the 
amount  available  per  individual  is 


NI  =  N 


_ bo _ 

I  ds  b(s ) 


where  &o  corresponds  to  the  biomass  of  a  single  individual  and  the  notation  (b)  =  f  ds  b 
has  been  used.  Growth  will  occur  when  Nj  >  Ncra  but  saturates  to  a  maximum  value  gm 
when  the  resources  are  plentiful.  For  example, 


9  —  9m 


A Ti  -  Ncrit 
Nj 


for  Nj  >  Ncrit, 


which  results  in  a  logistic-like  system: 


9m 


,  d2 

b  +  m,-—r;b 
osz 


(3) 
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The  coefficients  gm  (maximum  growth  rate),  d  (death  rate),  and  bc  =  Nbo/Ncrit  (carrying 
capacity)  are  functions  of  s.  For  this  form, 

K(s\E)  =  K(s\(b))  =  gm  (f  -  ■£)  -  d. 

The  environment  here  is  the  density  of  competitors  for  the  resource  as  well  as  the  value  of 
N  implicit  in  the  carrying  capacity  bc. 


3  Evolutionarily  Stable  Strategies  (ESS) 


If  the  system  modeled  by  Equation  3  is  analyzed  without  considering  mutation  (m  =  0), 
there  are  many  possible  singular  solutions  of  the  form 


b(s,  t )  =  b(t)S(s  —  s ) 


with  b  satisfying 


p  =  K(s\b)b 


since  E  =  (b)  =  b.  This  results  in  just  the  logistic  equation  and  steady  solutions  where 
TZ(s\b)  =  0  indicating 

d{s) 


b  bc(s)  1  - 


9m.{s , 


Let  us  introduce  a  different  organism  with  phenotype  s':  b  =  b(t)S(s  —  s)  +  b'(t)S(s  —  s'); 
if  its  biomass  is  very  low  when  compared  to  b.  Equation  3  becomes 


-  nm  = 


9m(s')  (  1  - 


bc(s') 


-  d(s') 


and  when  the  coefficients  are  time- independent,  we  can  write  this  in  the  simpler  form 


ld_ 

br  dt 


b' 


9m(s') 


b(s')  -  b(s) 

bc(s') 


where  b(sr)  is  the  equilibrium  population  for  s  =  s' .  For  this  model,  the  population  with  the 
highest  equilibrium  value  will  exclude  all  others:  b'  will  decay  if  b(s')  <  b(s).  Alternatively, 
the  ESS  maximizes  bc(s)[gm(s )  —  d(s)]/gm(s),  and  it’s  also  the  one  which  can  survive  on 
the  minimum  resource  Nj. 

In  other  situations,  there  may  be  no  single  species  which  can  out-compete  all  others. 
For  example,  consider  a  case  in  which  different  species  use  somewhat  different  resources  so 
that  Nj  for  a  species  s  depends  on  the  nearby  biomass  but  not  species  far  away  in  trait 
space.  E.g,  is  s  represents  the  size  of  a  herbivore,  it  will  generally  feed  on  a  limited  range 
of  plant  sizes.  This  kind  of  local  competition  could  look  like 


E  —  gm 


ds'w(s\s')b(s') 


d 


with  w  a  peaked  function.  Since  we  will  take  bc  to  be  constant,  we  can  just  choose  it  to  be 
one  hereafter;  w  will  be  normalized  to  have  unit  integral.  Figure  1  shows  an  example  with 
speciation  events  and  a  final  state  with  6  “species.” 


136 


1 


0.8 


0.6 


0.4 


0.2 


0 

0  50  100  150  200  250 

Figure  1:  Biomass  density  vs.  time  for  the  case  of  local  competition.  The  y  axis  is  s;  the 
x  axis  is  time,  but  on  an  uneven  scale.  The  black  lines  show  where  the  time  scale  changes 
from  1  day  per  unit  to  5,  10,  50,  and  500  days. 

4  Adaptive  dynamics 

Dieckmann  points  out  that  the  dynamics  may  be  approximated  by  equations  for  the  net 
population  and  the  rate  of  movement  in  trait  space. 

We  can  look  at  the  (local)  mean  biomass  and  the  mean  value  of  s 

( b )  =  fsodsb(s ) 

<->  =  (sb)/(b) 

Applying  the  local  average  to  a  system  defined  in  Equation  2  and  ignoring  fluxes  out  of  the 
range  by  mutation,  we  get 

^L(b)  =  (n(s,b)b) 

which  is 

-K{(s),b)(b)  +  ^K"a2(b) 

from  Taylor  expansion.  For  evaluating  the  competition/  environmental  effects,  we  can  treat 
the  biomass  distribution  as  a  set  of  singular  values 

b  =  ( b)i5(s  -  (s)i) 
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Taking  the  local  s  moment  and  expanding  gives 


k{s>(6)  =  (sKb) 

=  ((» -  <‘)m 

~  lZ'cr2(b) 


The  speed  at  which  the  peaks  migrate  towards  the  optimal  value  depends  on  the  gradient 
and  on  the  width  of  the  distribution  a.  Spread  arises  from  mutation  or  from  migration 
mixing.  We  could  derive  an  equation  for  <r;  it  tends  to  increase  because  of  mn  and  decrease 
because  of  negative  curvature  in  7Z  which  acts  to  sharpen  peaks. 

As  a  simple,  and  analytically  tractable  example  for  a  single  peak,  consider  a  system 
with 


-R2(s  -  s0) 


2 


(b) 


where  i?o,  R2,  and  so  could  be  functions  of  time  as  the  population  or  other  processes  alter 
the  resource.  We  can  find  a  time-dependent  solution  (without  needing  to  approximate) 


b  =  bo(t)  exp 


1  [s  —  x(t )]2 

2  u2(f) 


when 


/  d 

rci—a  = 
dt 


1„*  q  m 

-R*crs  H - 

2  a 


d_ 

dt 


X  =  —R2<JZ(x  —  So) 


1  d  _  1  2  m 

--b0-R0--R2(x-s0)  -  — 


with  Rj  =  Rj(  1  —  crfeo)- 

The  center  of  the  distribution  moves  towards  the  peak  in  the  growth  rate  curve  at  a  speed 
which  depends  of  the  slope  at  the  current  center  location  [— R%{x  —  so)]  and  the  variance 
of  the  distribution.  The  width  asymptotes  to  the  value  {2m/ i?^)1/4  (if  1Z"  is  negative) 
which  narrows  as  the  growth  rate  curve  becomes  more  sharply  peaked.  The  population  will 
increase  or  decrease  depending  on  the  sign  of  Rq  —  m/a2',  the  dynamics  of  the  resource  will 
adjust  Rq  until  it  reaches  the  equilibrium  value. 

The  environment  will  be  time-dependent  because  of  many  external  factors,  so  that  1Z 
itself  has  a  whole  spectrum  of  variability.  The  velocity  for  movement  in  trait  space  will 
therefore  be  fluctuating,  and  we  cannot  expect  the  system  to  be  in  equilibrium.  If  the 
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external  changes  are  very  slow  compared  to  mutation  times,  the  system  will  be  very  close 
to  the  ESS  for  current  conditions: 


dlZ(s,  t ) 


dsj 


s=s(t) 


~  0 


In  contrast,  if  the  environmental  fluctuations  are  fast,  the  center  of  the  distribution  will  not 
be  able  to  keep  up,  and  the  mean  trait  value  will  settle  at  the  point  where  the  time-averaged 
velocity  is  zero 


dlZ(s,  t ) 


dsj 


~  0 


For  intermediate  time-scales,  comparable  to  mutation  times,  the  center  of  the  distribution 
will  partially,  but  not  completely,  track  the  variations.  The  population  dynamics  is  implicit 
here,  basically  maintaining  as  we  shall  see,  it  may  play  a  much  more  active  role 

with  multiple  trophic  levels. 


5  Time-dependent 

We  can  also  consider  a  system  distributed  in  trait-  and  physical  space  in  which  a  changing 
environment  produces  time-dependent  forcing.  A  simple  way  to  do  this  is  to  change  the 
reproduction  rate  in  such  a  way  that  favors  different  phenotypes  at  different  times.  This 
could  model  such  changing  environmental  parameters  as  mixed  layer  depth,  or  sunlight, 
for  example.  During  the  portion  of  the  forcing  cycle  when  a  certain  phenotype  is  not  fit, 
mutation  can  replenish  such  populations  which  would  otherwise  slowly  die  out.  We  can 
then  see  a  seasonal  cycle  in  the  community  structure  as  well  as  the  biomass.  Alternatively, 
selective  grazing  in  which  the  zooplankton  feed  more  heavily  on  the  more  abundant  species 
can  also  prevent  the  non-optimal  species  from  disappearing. 

If  the  growth  rate  is  time-dependent,  ^  0,  a  time-average  of  the  biomass  can  be 
taken  to  find 

f  ln17W)  =  =  R(S'\^) 

Then,  based  on  this  result,  s  will  represent  an  evolutionarily  stable  strategy  if 

i?(s|s)  =  0  [definition]  and  R(s'\s)  <0  for  s'  /  s 

which  indicates  that  organisms  with  a  non-optimal  phenotype  will  die  out  representing 
natural  selection.  Locally,  these  conditions  become 

n  n2 

-^R(s'\s)\s,=.  =  0  ,  -^R(s'\s)\s,=-<0, 

but  since  the  competition  is  temporally  local,  it  is  possible  that  several  species  could  coexist 
with  one  or  another  dominant  at  different  times. 


139 


0.8 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

0 

0  50  100  150  200  250  300  350 

t 

Figure  2:  Biomass  density  vs.  time  for  s  =  0.41  (late)  and  for  the  case  of  local  competition. 
The  y  axis  is  s;  the  x  axis  is  time,  but  on  an  uneven  scale.  The  black  lines  show  where  the 
time  scale  changes  from  1  day  per  unit  to  5,  10,  50,  and  500  days. 
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6  Relationship  to  stability  theory 

To  understand  the  relationship  to  stability  theory,  look  again  at  the  linearized  model 


with  the  basic  state  having 


|6'  =  K(,|5)6'  +  5||W 


b(s)'R.(s\b(s))  =  0 


Standard  problems:  b(s)  /  0  =>■  7 Z(s,b(s))  =  0,  the  first  term  vanishes,  and  we  deal 
with  the  second  term,  often  in  the  form 


d  ,  _  r  dnt 

dtbi  - bi  db, 


6' 


b=b 


Evolutionarily  Stable  Strategies:  b(s)  =  0  for  s  /  s  so  that  lZ(s\b)  ^  0  and  the  first 
term  is  the  important  one. 


In  the  case  with  mutation 


b1Z(s\b)  +  \72sm.b  =  0 


^b'  =  n{s\b)b'  +  b^(b')  +  V2smb' 
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both  terms  enter.  Near  the  “hot-spots”  where  b  is  large,  1Z  is  relatively  small,  the  second 
term  dominates  and  the  problem  looks  like  a  standard  stability  system  modified  by  diffusion. 
Far  away,  however,  b  is  small,  and  the  basic  state  has  a  decaying  form  with  1Z  nonzero.  If 
1Z  changes  and  becomes  positive  in  that  region,  the  population  can  “tunnel”  into  the  local 
maximum  and  grow.  The  amount  of  time  this  will  take  depends  on  distance  and  the 
mutation  rate. 


7  Physical  mixing 

We  now  consider  a  model  with  a  single  spacial  dimension,  y,  in  addition  to  the  single  phe¬ 
notype  dimension,  s.  Our  growth  rate  varies  in  physical  space,  and  the  biomass  also  diffuses 
in  physical  space,  but  we  neglect  the  diffusion  in  trait  space  that  represents  mutation. 

§ib  =  bTZ(s,y\{b))  +  K-^b, 

with 

n=[Ro-  \R2(s  -  S0(y))2]  [1  -  <&>]  -  d. 

Figure  3  shows  the  biomass  concentrates  into  a  limited  number  of  species  occupying  over¬ 
lapping  latitude  ranges;  intermediate  species  are  excluded. 


b 


S 


Figure  3:  Biomass  density  b(s)  vs.  y  (non-dimensional)  with  physical  mixing. 

Because  there  is  no  mutation,  we  can  analyze  this  behavior  by  examining  singular  solu¬ 
tions  with  the  single-phenotype  ansatz, 

b(s,y,t)  =  b(y,  t)5(s  -  s). 
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Applying  this  ansatz  yields  the  PDE  in  y  and  t  that  governs  6, 


d_ 

dt 


b=([R0- 


\R2{s  -  «o(y))2]  [!  -  b]  -d)b  + 


K 


Note  that  the  linear  part  of  this  operator  is  a  Schrodinger  equation, 


d 2 

iy)b  + 


5, 


which  has  a  trapped  mode  if  the  region  where  lZ(s,y  |0)  >  0  is  big  enough  and  strong  enough 
to  overcome  diffusive  losses  into  the  region  where  the  death  rate  dominates  (Kierstead  and 
Slobodkin). 

To  determine  whether  a  steady  single-phenotype  solution  is  linearly  stable  to  nearby 
species  (i.e.  is  a  local  ESS),  consider  small  perturbations,  b(s,y,t)  =  b(y)  +  b'(s  +  s',y,t ). 
To  linear  order  in  s'  and  b' ,  the  perturbations  are  governed  by 


_  <97?  _  <92 

b'K{s,  y,  b )  +  s'b—(s,  y,  b)  +  ■ 


If  s  is  a  local  ESS,  the  growth  rate  for  b'  will  be  zero,  and  we  would  need  one  higher  order 
is  s'  to  decide  if  the  state  is  a  local  minimum  or  maximum.  Thus,  we  want  to  find  the  s 
such  that 

_  <97?  _  <9 2 

0  =  b'U(s,  y ,  b )  +  s'b—(s,  y.  b)  +  K-^b' . 

in  the  limit  as  s'  —>  0.  Multiplying  this  equation  by  b  and  the  equation  for  the  latter  by  b' 
and  integrating  over  space  implies 


,  -  , ,  dlZ  .  — . 

dx  bb  — ( s,y,b )  =  0 


In  the  limit,  b'  — >  b ,  so  we  can  pick  an  s,  hnd  b,  evaluate  this  integral,  and  search  for  the 
s  where  it  vanishes.  Species  near  this  value  will  be  excluded;  however,  the  growth  rate  for 
ones  far  away,  found  from 

d2 

ab'  =  TZ(s ,  y,  b)b'  +  n— ^b' , 

may  be  positive.  In  that  case,  another  band  will  appear.  This  can  happen  because  b  varies 
spatially,  so  it  is  not  effective  at  excluding  others  in  the  regions  where  b  is  small  and  is  dying 
out  but  sustained  by  diffusion  from  the  positive  growth  regions. 

Simulations  of  this  spatial-diffusion  model  give  the  following  qualitative  results: 

•  The  single  species  solution  has  reduced  amplitude  and  will  die  out  if  k  is  big  and  the 
volume  average  of  72.(s|0)  is  negative. 

•  If  b(s,  y)  is  large  enough,  single  solutions  may  damp  the  rest  of  the  held,  even  when 
it  would  otherwise  grow. 

•  Discrete  species  with  different  phenotypes  emerge  in  spatially-separated  bands. 
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8  Predator-prey 


We  now  present  a  model  with  a  predator  species  (zooplankton,  Z)  and  a  prey  species  (phy¬ 
toplankton,  P)  evolving  in  a  one-dimensional  trait  space.  The  prey  has  a  trait-dependent 
birth  rate,  the  predator  has  a  constant  death  rate,  and  the  predator  feeds  on  the  prey 
according  to  their  respective  traits,  as  governed  by  the  interaction  kernel,  G. 


D 

ITt 

D 

m 


P(s)  = 
Z(s')  = 


n(s)N  -  /  ds'  G(s  -  s')Z(s') 


P 


a  /  ds  G(s  —  s')P(s)  —  dz 


Z, 


with  /x  peaked  at  s 
be  preserved: 


0  and  G(s  —  s')  peaked  at  zero.  We  also  require  that  the  total  biomass 


N+  ds  P(s)  +  /  ds'  Z(s'). 


This  system  has  a  singular  solution, 


(PZ)  =  (P5(s),Z5(s')), 


which  is  an  ESS,  but  it  may  be  unstable  if 


|G"(0)|Z>|//(0)|iV. 


or  if 

G(0)  MO)  ’ 

so  that  the  predation  function  is  sharply  curved  and  inhibits  growth  only  in  a  small  neigh¬ 
borhood  of  the  single-phenotype  solutions.  The  figures  show  two  examples  of  P(s,  t )  and 
Z(s't)  for  weakly  and  more  strongly  unstable  situations. 


9  Advection 

Another  way  to  make  the  environment  change  periodically  is  to  make  the  time  derivative 
into  an  advective  derivative,  and  to  impose  a  flow  that  models  ocean  circulation.  In  the 
model  that  was  simulated,  the  “ocean”  is  wind-driven,  two  layers  with  warm  water  above 
the  thermocline  and  cold  water  below),  and  quasi-geostrophic.  Figure  6  shows  a  snapshot 
of  the  mean  s  value  from  such  a  model  with  a  single  nutrient  and  a  single  predator.  The 
diversity  is  highest  in  the  region  where  the  western  boundary  currents  meet  and  proceed 
offshore  as  an  eddying  jet. 
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Figure  5:  Predator-prey  biomass  densities  for  G" /G  =  n"  /  ji. 
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Figure  6:  Mean  value  of  s  in  a  two-gyre,  eddying  model. 
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